home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Technotools
/
Technotools (Chestnut CD-ROM)(1993).ISO
/
lang_oth
/
complex
/
complex.lib
Wrap
Text File
|
1987-11-08
|
5KB
|
167 lines
{***********************************************************************}
{ }
{ TURBO Pascal library of Complex Number routines adapted from : }
{ October 1984 Dr. Dobbs Journal by John Lucas 30 Sept. 1984 }
{ }
{-----------------------------------------------------------------------}
{ }
{ Global Declarations needed by this Library : }
{ }
{ type }
{ complex = record }
{ re : real; }
{ im : real; }
{ end; }
{ }
{***********************************************************************}
procedure cadd (var result : complex; arg1,arg2 : complex);
begin
result.re := arg1.re + arg2.re;
result.im := arg1.im + arg2.im
end;
procedure csub (var result : complex; arg1,arg2 : complex);
begin
result.re := arg1.re - arg2.re;
result.im := arg1.im - arg2.im
end;
procedure cmult (var result : complex; arg1,arg2 : complex);
begin
result.re := arg1.re * arg2.re - arg1.im * arg2.im;
result.im := arg1.im * arg2.re - arg1.re * arg2.im
end;
procedure cdiv (var result : complex; arg1,arg2 : complex);
var
denom : real;
begin
denom := sqr(arg2.re) + sqr(arg2.im);
result.re := (arg1.re * arg2.re + arg1.im * arg2.im)/denom;
result.im := (arg1.im * arg2.re - arg2.re * arg2.im)/denom
end;
procedure polar (arg : complex; var modulus,amplitude : real);
const
lnmaxreal = 87.49823353;
piover2 = 1.570796327;
closest = 1E-19;
begin
with arg do
begin
if abs(re) < closest then
re := 0.0;
if abs(im) < closest then
im := 0.0;
modulus := sqrt(sqr(re) + sqr(im));
if im = 0.0 then
amplitude := 0.0
else
if re = 0.0 then
if im = 0.0 then
amplitude := piover2
else
amplitude := -piover2
else
if (ln(abs(im)) - ln(abs(re)) > lnmaxreal) then
if re > 0.0 then
if im > 0.0 then
amplitude := piover2
else
amplitude := -piover2
else
if im > 0.0 then
amplitude := -piover2
else
amplitude := piover2
else
amplitude := arctan(im/re)
end;
end;
procedure ctopower (var result : complex; arg : complex; power : integer);
var
i : integer;
modulus,amplitude,newmod,powamp : real;
begin
if power = 0 then
begin
result.re := 1.0;
result.im := 0.0
end
else
begin
polar(arg,modulus,amplitude);
newmod := 1.0;
if power > 0 then
for i := 1 to power do
newmod := newmod * modulus
else
for i := 1 to abs(power) do
newmod := newmod/modulus;
powamp := power * amplitude;
result.re := newmod * cos(powamp);
result.im := newmod * sin(powamp)
end;
end;
procedure cexp (var result : complex; arg : complex);
var
expo : real;
begin
expo := exp(arg.re);
result.re := expo * cos(arg.im);
result.im := expo * sin(arg.im)
end;
procedure cln (var result : complex; arg : complex);
var
modulus,amplitude : real;
begin
polar(arg,modulus,amplitude);
result.re := ln(modulus);
result.im := amplitude
end;
procedure ctoc (var result : complex; arg1,arg2 : complex);
var
logpart,expo : complex;
begin
cln(logpart,arg1);
cmult(expo,arg2,logpart);
cexp(result,expo)
end;
procedure csin (var result : complex; arg : complex);
var
exp1,exp2,part1,part2,sum,divisor : complex;
begin
exp1.re := -arg.im;
exp1.im := arg.re;
cexp(part1,exp1);
exp2.re := arg.im;
exp2.im := -arg.re;
cexp(part2,exp2);
csub(sum,part1,part2);
divisor.re := 0.0;
divisor.im := 2.0;
cdiv(result,sum,divisor)
end;
procedure ccos (var result : complex; arg : complex);
var
exp1,exp2,part1,part2,sum,divisor : complex;
begin
exp1.re := -arg.im;
exp1.im := arg.re;
cexp(part1,exp1);
exp2.re := arg.im;
exp2.im := -arg.re;
cexp(part2,exp2);
cadd(sum,part1,part2);
divisor.re := 2.0;
divisor.im := 0.0;
cdiv(result,sum,divisor)
end;